Regional Poverty Convergence across Districts in Indonesia:A Distribution Dynamics Approach

Suggested citation: > Miranti, Ragdad Cani.(2020). Regional Poverty Convergence across Districts in Indonesia: A Distribution Dynamics Approach https://rpubs.com/canimiranti/distribution_dynamics_poverty514districts

This work is licensed under the Creative Commons Attribution-Non Commercial-Share Alike 4.0 International License.

1 Original data source

Data of Headcount Index (Poverty Rate) are derived from the Indonesia Central Bureau of Statistics (Badan Pusat Statistik Republik of Indonesia). https://www.bps.go.id/

Acknowledgment:

Material adapted from multiple sources, in particular from Magrini (2007).

2 Libraries

knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(skimr)
library(kableExtra)    # html tables 
library(pdfCluster)    # density based clusters
library(hdrcde)        # conditional density estimation 
library(plotly)

library(intoo)
library(barsurf)
library(bivariate)

library(np)
library(quantreg)

library(basetheme)
basetheme("minimal")

library(viridis)
library(ggpointdensity)
library(isoband)

#library(MASS)
library(KernSmooth)


# Change the presentation of decimal numbers to 4 and avoid scientific notation
options(prompt="R> ", digits=4, scipen=7)

3 Tutorial objectives

• Study the dynamics of univariate densities

• Compute the bandwidth of a density

• Study mobility plots

• Study bi-variate densities

• Study density-based clustering methods

• Study conditional bi-variate densities

4 Import data

dat <- read_csv("poverty.csv")
dat <- as.data.frame(dat)

5 Transform data

Since the data is in relative terms, let us rename the variables and add new variables.

dat <- dat %>% 
  mutate(
    rel_pov2010 = pov2010/mean(pov2010),
    rel_pov2018 = pov2010/mean(pov2018)
  )
dat

6 Descriptive statistics

skim(dat)

── Data Summary ────────────────────────
                           Values
Name                       dat   
Number of rows             514   
Number of columns          13    
_______________________          
Column type frequency:           
  character                1     
  numeric                  12    
________________________         
Group variables            None  

── Variable type: character ────────────────────────────────────────────────────────────
  skim_variable n_missing complete_rate   min   max empty n_unique whitespace
1 District              0             1     4    26     0      514          0

── Variable type: numeric ──────────────────────────────────────────────────────────────
   skim_variable n_missing complete_rate    mean       sd      p0      p25      p50
 1 Code                  0             1 4558.   2678.    911     1803.    3550    
 2 pov2010               0             1   15.6     9.42    1.67     9.01    13.3  
 3 pov2011               0             1   14.6     8.92    1.5      8.14    12.4  
 4 pov2012               0             1   13.9     8.55    1.33     7.79    11.8  
 5 pov2013               0             1   13.8     8.59    1.75     7.76    11.7  
 6 pov2014               0             1   13.1     7.93    1.68     7.38    11.2  
 7 pov2015               0             1   13.5     8.25    1.69     7.59    11.6  
 8 pov2016               0             1   13.2     8.20    1.67     7.40    11.2  
 9 pov2017               0             1   13.0     7.98    1.76     7.35    11.1  
10 pov2018               0             1   12.3     7.84    1.68     6.99    10.2  
11 rel_pov2010           0             1    1       0.603   0.107    0.577    0.854
12 rel_pov2018           0             1    1.27    0.763   0.135    0.730    1.08 
       p75    p100 hist 
 1 7174.   9471    ▇▆▂▇▃
 2   19.6    49.6  ▇▇▂▁▁
 3   18.9    47.4  ▇▇▂▁▁
 4   17.8    45.9  ▇▇▂▁▁
 5   17.4    47.5  ▇▆▂▁▁
 6   16.6    44.5  ▇▆▂▁▁
 7   16.9    45.7  ▇▇▂▁▁
 8   16.3    45.1  ▇▇▂▁▁
 9   16.0    43.6  ▇▇▂▁▁
10   15.1    43.5  ▇▆▂▁▁
11    1.25    3.18 ▇▇▂▁▁
12    1.59    4.02 ▇▇▂▁▁

xy <- dat %>% 
select(
pov2010,
pov2018,
  ) %>% 
  mutate(
    x = pov2010,
    y = pov2018
  ) %>% 
  select(
    x,
    y
  )

7 Univariate dynamics

7.1 Select bandwiths

Select bandwidth based on function dpik from the package KernSmooth

h_rel_pov2010<- dpik(dat$rel_pov2010)
h_rel_pov2010

[1] 0.1026

h_rel_pov2018 <- dpik(dat$rel_pov2018)
h_rel_pov2018

[1] 0.1299

7.2 Plot each density

dis_rel_pov2010 <- bkde(dat$rel_pov2010, bandwidth = h_rel_pov2010)
dis_rel_pov2010 <- as.data.frame(dis_rel_pov2010)
ggplot(dis_rel_pov2010, aes(x, y)) + geom_line() + 
  theme_minimal() 

dis_rel_pov2018 <- bkde(dat$rel_pov2018, bandwidth = h_rel_pov2018)
dis_rel_pov2018 <- as.data.frame(dis_rel_pov2018)
ggplot(dis_rel_pov2018, aes(x, y)) + geom_line() + 
  theme_minimal() 

7.3 Plot both densities

There are two methods for plot both densities,i.e. Kernsmooth and bandwith default of ggplot. I prefer using the gglot ( Method 2).

8 Method 2 (ggplot)

Using the bandwidth default of ggplot Manual labels are not yet implemented in the ggplotly function

rel_pov2010 <- dat %>% 
  select(rel_pov2010) %>% 
  rename(rel_var = rel_pov2010) %>% 
  mutate(year = 2010)

rel_pov2018 <- dat %>% 
  select(rel_pov2018) %>% 
  rename(rel_var = rel_pov2018) %>% 
  mutate(year = 2018)

rel_pov2010pov2018 <- bind_rows(rel_pov2010, rel_pov2018)

rel_pov2010pov2018 <- rel_pov2010pov2018 %>% 
  mutate(year = as.factor(year))
head(rel_pov2010pov2018)

dis_rel_pov2010pov2018 <- ggplot(rel_pov2010pov2018,aes(x=rel_var, color=year)) +
  geom_density() + 
  theme_minimal() 
dis_rel_pov2010pov2018

Using plotly

ggplotly(dis_rel_pov2010pov2018)

9 Bivariate density

9.1 Using Mobility scatterplot

dat %>% 
  ggplot(aes(x = rel_pov2010, y = rel_pov2018)) +
  geom_point(alpha=0.5) +
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept = 1, linetype="dashed") + 
  geom_vline(xintercept = 1, linetype="dashed") +
  theme_minimal() +
  labs(subtitle = "Relative Pov2018",
       x = "Relative Pov2010",
       y = "") +
  theme(text=element_text(family="Palatino")) 

Fit a non-linear function

dat %>% 
  ggplot(aes(x = rel_pov2010, y = rel_pov2018)) +
  geom_point(alpha=0.5) + 
  geom_smooth() + 
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept = 1, linetype="dashed") + 
  geom_vline(xintercept = 1, linetype="dashed") +
  theme_minimal() +
  labs(subtitle = "Relative pov2018",
       x = "Relative pov2010",
       y = "") +
  theme(text=element_text(family="Palatino")) 

Not that the nonlinear fit crosses the 45-degree line two times from above.

9.2 Using the Bivariate package

bivariate <- kbvpdf(dat$rel_pov2010, dat$rel_pov2018, h_rel_pov2010, h_rel_pov2018) 

plot(bivariate,
      xlab="Relative Poverty 2010", 
      ylab="Relative Poverty 2018")
abline(a=0, b=1)

abline(h=1, v=1)

plot(bivariate,
      TRUE,
      xlab="Relative Poverty 2010", 
      ylab="Relative Poverty 2018")

9.3 Using ggplot (stat_density_2d)

dat %>% 
  ggplot(aes(x = rel_pov2010, y = rel_pov2018)) +
  geom_point(color = "lightgray") + 
  geom_smooth() + 
  #geom_smooth(method=lm, se=FALSE) + 
  stat_density_2d() +
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept = 1, linetype="dashed") + 
  geom_vline(xintercept = 1, linetype="dashed") +
  theme_minimal() +
  labs(subtitle = "Relative pov2018",
       x = "Relative pov2010",
       y = "") +
  theme(text=element_text(family="Palatino")) 

dat %>% 
  ggplot(aes(x = rel_pov2010, y = rel_pov2018)) +
        stat_density_2d(aes(fill = stat(nlevel)), geom = "polygon") + 
  scale_fill_viridis_c() +
        geom_abline(aes(intercept = 0, slope = 1)) +
        geom_hline(yintercept = 1, linetype="dashed") + 
        geom_vline(xintercept = 1, linetype="dashed") + 
  theme_minimal() +
        labs(x = "Relative Poverty Rate 2010",
             y = "Relative Poverty Rate 2018") +
        theme(text=element_text(size=8, family="Palatino"))

9.4 Using Stochastic Kernel Package

There are two ways of analyzing the use of stochastic kernel package: 1. Contour Plots and 2. Surface Plots

9.4.1 Contour Plots

pov2010pov2018 <- cbind(dat$rel_pov2010, dat$rel_pov2018)
pov2010pov2018_dis <- bkde2D(pov2010pov2018, bandwidth = c(h_rel_pov2010, h_rel_pov2018)) 

contour(pov2010pov2018_dis$x1,pov2010pov2018_dis$x2,pov2010pov2018_dis$fhat)
abline(a=0, b=1)

9.4.2 Surface Plots

plot_ly(x=pov2010pov2018_dis$x1, y=pov2010pov2018_dis$x2, z=pov2010pov2018_dis$fhat) %>% add_surface()

plot_ly(x=pov2010pov2018_dis$x1, y=pov2010pov2018_dis$x2, z=pov2010pov2018_dis$fhat, type = "contour", contours = list(showlabels = TRUE))  %>%
  colorbar(title = "Density")

10 Conditional density analysis

10.1 Using the hdrcde package

pov2010pov2018_cde <- cde(dat$pov2010, dat$pov2018)

Increase the number of intervals to 20

pov2010pov2018_cde2 <- cde(dat$pov2010, dat$pov2018, nxmargin = 20)

plot(pov2010pov2018_cde)

plot(pov2010pov2018_cde2)

High density regions

plot(pov2010pov2018_cde, plot.fn="hdr")
abline(a=0, b=1)

abline(h=1, v=1)

11 References

• Magrini, S. (2007). Analysing convergence through the distribution dynamics approach: why and how?. University Ca’Foscari of Venice, Dept. of Economics Research Paper Series No, 13.

• Mendez C. (2020). Classical sigma and beta convergence analysis in R: Using the REAT 2.1 Package. R Studio/RPubs. Available at https://rpubs.com/quarcs-lab/classical-convergence-reat21

• Mendez C. (2020). Univariate distribution dynamics in R: Using the ggridges package. R Studio/RPubs. Available at https://rpubs.com/quarcs-lab/univariate-distribution-dynamics

• Mendez, C. (2020) Regional efficiency convergence and efficiency clusters. Asia-Pacific Journal of Regional Science, 1-21.

• Mendez, C. (2019). Lack of Global Convergence and the Formation of Multiple Welfare Clubs across Countries: An Unsupervised Machine Learning Approach. Economies, 7(3), 74.

• Mendez, C. (2019). Overall efficiency, pure technical efficiency, and scale efficiency across provinces in Indonesia 1990 and 2010. R Studio/RPubs. Available at https://rpubs.com/quarcs-lab/efficiency-clusters-indonesia-1990-2010

• Mendez-Guerra, C. (2018). On the distribution dynamics of human development: Evidence from the metropolitan regions of Bolivia’’. Economics Bulletin, 38(4), 2467-2475.

END